get_survival_curve <- function(var, nice_name) {
  
  mod <- coxph(as.formula(paste0("Surv(tdfs, deltadfs) ~ ", var)),
     data = bmt)
  
  if (is.factor(bmt[[var]])) {
    newdat <- tibble(
    !!var := unique(bmt[[var]])
  )
  }
  else {
    newdat <- tibble(
    !!var := round(quantile(bmt[[var]],probs= c(.25,.75), na.rm=TRUE),2))
    
  }
  
  
  est_by_group <- pmap_df(newdat, ~{
    ggsurvfit(survfit(mod, newdata =tibble(!!var := .x)))$data %>%
      mutate(strata = .x)}
    )
  
  if(!is.factor(bmt[[var]]))  {
    est_by_group <- est_by_group %>%
      mutate(strata=paste0(nice_name, " = ", strata))
  }

  est_by_group %>%
    select(estimate, strata, time, conf.high, conf.low) %>%
    dplyr::arrange(strata, time) |>
    dplyr::group_by(strata) |>
    dplyr::mutate(
    xmin = time,
    xmax = dplyr::lead(time)   # next step’s start
  ) |>
    mutate(strata=gsub(paste0(var,"="), 
                       "",
                       strata)) %>%
    ggplot(aes(x=time,y=estimate,color=strata)) +
    geom_step() +
    geom_rect(aes(xmin = xmin, xmax = xmax, 
                ymin = conf.low, ymax = conf.high,
                fill=strata),
              alpha = 0.15,
              inherit.aes=FALSE) +
    theme_c(legend.text=element_text(size=16),
            plot.title=element_text(size=16,
                                    margin = margin(b = 5))) +
    viridis::scale_color_viridis(discrete=TRUE,
                                 option="magma", end=.8, begin=.3) +
    viridis::scale_fill_viridis(discrete=TRUE,
                                 option="magma", end=.8, begin=.3) +
    labs(y="Disease-Free Survival", x= "Time (Days)",
         title=paste0("By ", nice_name), color="", fill="") +
    scale_x_continuous(n.breaks=7) +
    guides(color=guide_legend(override.aes = list(linewidth=2)))
  
}

pltlist <- map2(covariates,nice_names,get_survival_curve)

# package for arranging plots
library(patchwork)

# arrange plots into a grid 
(pltlist[[1]] + pltlist[[2]])/ (pltlist[[3]] + pltlist[[4]]) / (pltlist[[5]] + pltlist[[6]]) / (pltlist[[7]] + pltlist[[8]])/(pltlist[[9]] + pltlist[[10]]) + plot_annotation(
  title = "Survival Curves by Baseline Variables of Interest",
  theme=theme_c(plot.title=element_text(face="bold",hjust=.5, size=20)))

# save plot 
ggsave(here("figures/surv_by_baseline.jpeg"))

library(gt)
library(gtsummary)

get_model_summary <- function(var, nice_name, model) {
  
  # model <- coxph(as.formula(paste0("Surv(tdfs, deltadfs) ~ ", var)),
  #    data = bmt)
  
  if (var == "male" | var == "donorcmv" | var == "donormale" ) {
    model_summary <- tbl_regression(model, exponentiate=TRUE,
                                  conf.level=.9,
                 estimate_fun=label_style_sigfig(digits=3),
                 pvalue_fun = label_style_pvalue(digits = 3),
                 show_single_row=all_of(var),
                 label=setNames(
                   list(nice_name, paste0("tt(", nice_name, ")")), 
                   c(var, paste0("tt(",var, ")"))))
  } else {
     model_summary <- tbl_regression(model, exponentiate=TRUE,
                                  conf.level=.9,
                 estimate_fun=label_style_sigfig(digits=3),
                 pvalue_fun = label_style_pvalue(digits = 3),
                 label=setNames(
                   list(nice_name, paste0("tt(", nice_name, ")")), 
                   c(var, paste0("tt(",var, ")"))))
  }
  
 
  return(model_summary)
  
}

get_model <- function(var) {
  
  model <- coxph(as.formula(paste0("Surv(tdfs, deltadfs) ~ ", var)),
     data = bmt)
  
  return(model)
  
}

model_list <- map(covariates, get_model)
names(model_list) <- covariates


# linear spline basis 1: min(t, c) # basis 2: (t - c)+)
#----------------------------
# plot Schoenfeld residuals 
#----------------------------

library(splines)

get_res_plot <- function(var, nice_name, model) {
  
  res <- as.data.frame(cox.zph(model, transform="km")$y)[[var]]
  # pval <- cox.zph(model, transform="km")$table[var,"p"] %>%
  #   signif(3)
  # 
  if(is.null(res)) {return(NULL)}
  time <- cox.zph(model, transform="km")$time
  df <- tibble(var = var,
         time = time,
          y=res) 
  # glimpse(df)
  
  spline_fit <- lm(y ~ ns(time, df = 3), data = df)
  ci <- predict(spline_fit, newdata = df,
                interval = "confidence",
                level = 0.95)
  df_ci <- bind_cols(df, as.data.frame(ci)) 
  
  df_ci %>%
    mutate(
      spline = predict(spline_fit)
    ) %>%
      ggplot(aes(x=time)) +
       geom_point(aes(y=y),alpha=.3) +
       geom_line(aes(x=time,y=upr), linetype=2, linewidth=1.05) +
       geom_line(aes(x=time,y=lwr), linetype=2, linewidth=1.05) +
       geom_line(aes(x=time,y=spline),linewidth=1.05) +
    theme_c(plot.title=element_text(face="plain", size=10,hjust=0)) +
    labs(title=paste0(
                      "\nSchoenfeld Residuals for: ", nice_name,
                     # "\nP-value: ", pval,
                       "\n"))
    
}

plotlist <- pmap(list(var=covariates,
          nice_name=nice_names,
          model=model_list), 
     get_res_plot)


cowplot::plot_grid(plotlist=plotlist,ncol=2)

#-------------------------------------------------------------
# add time-varying covariates based on residual plots above 
#-------------------------------------------------------------

# linear interaction with time 
model_list[["donorage"]] <- coxph(
  Surv(tdfs, deltadfs) ~ donorage + tt(donorage),
     data = bmt,
      tt=function(x,t,...) x*t )


# piecewise linear interaction with time (to handle U-shapes)
# before the cutpoint, slope is beta_1
# after the cutpoint, slope is beta_2
# continuous at cutpoint 
model_list[["mtx"]] <- coxph(
  Surv(tdfs, deltadfs) ~ mtx + tt(mtx),
     data = bmt %>% mutate(mtx = ifelse(
       mtx=="Prophylactic Use of Methotrexate", 1, 0)),
      tt=function(x,t,...) cbind( x * pmin(t, 600), x * pmax(0, t - 600))) 

get_model_exclude_outlier <- function(var) {
  
  model <- coxph(as.formula(paste0("Surv(tdfs, deltadfs) ~ ", var)),
     data = bmt %>% filter(id != 69))
  
  return(model)
  
}

model_list_exclude_outlier <- map(covariates, get_model_exclude_outlier)
names(model_list_exclude_outlier) <- covariates


# linear spline basis 1: min(t, c) # basis 2: (t - c)+)
#----------------------------
# plot Schoenfeld residuals 
#----------------------------

library(splines)

plotlist <- pmap(list(var=covariates,
          nice_name=nice_names,
          model=model_list_exclude_outlier), 
     get_res_plot)


cowplot::plot_grid(plotlist=plotlist,ncol=2)

# MTX
mtx_mod <- coxph(
  Surv(tdfs, deltadfs) ~ mtx + tt(mtx),
     data = bmt %>% mutate(mtx = ifelse(
       mtx=="Prophylactic Use of Methotrexate", 1, 0)),
      tt=function(x,t,...) x*t)

tbl_regression(mtx_mod, exponentiate=TRUE, conf.level=.9)

# FAB
fab_mod <- coxph(
  Surv(tdfs, deltadfs) ~ fab + tt(fab),
     data = bmt %>% mutate(fab = ifelse(
       fab=="Other", 0, 1)),
      tt=function(x,t,...) x*t)

tbl_regression(fab_mod, exponentiate=TRUE, conf.level=.9)

# cmv
cmv_mod <- coxph(
  Surv(tdfs, deltadfs) ~ cmv + tt(cmv),
     data = bmt %>% mutate(cmv = ifelse(
       cmv=="CMV Positive", 1, 0)),
      tt=function(x,t,...) x*t)

tbl_regression(cmv_mod, exponentiate=TRUE, conf.level=.9)

# donor cmv
cmv_mod <- coxph(
  Surv(tdfs, deltadfs) ~ donorcmv + tt(donorcmv),
     data = bmt %>% mutate(donorcmv = ifelse(
       donorcmv=="CMV Positive", 1, 0)),
      tt=function(x,t,...) x*t)

tbl_regression(cmv_mod, exponentiate=TRUE, conf.level=.9)

# sex 
sex_mod <- coxph(
  Surv(tdfs, deltadfs) ~ male + tt(male),
     data = bmt %>% mutate(male = ifelse(
       male=="Male", 1, 0)),
      tt=function(x,t,...) x*t)

tbl_regression(sex_mod, exponentiate=TRUE, conf.level=.9)

# donor sex 
sex_mod <- coxph(
  Surv(tdfs, deltadfs) ~ donormale + tt(donormale),
     data = bmt %>% mutate(donormale = ifelse(
       male=="Male", 1, 0)),
      tt=function(x,t,...) x*t)

tbl_regression(sex_mod, exponentiate=TRUE, conf.level=.9)

# donor age 
donorsex_mod <- coxph(
  Surv(tdfs, deltadfs) ~ donorage + tt(donorage),
     data = bmt,
      tt=function(x,t,...) x*t)

tbl_regression(donorsex_mod, exponentiate=TRUE, conf.level=.9)

# disease group  
disgroup_mod <- coxph(
  Surv(tdfs, deltadfs) ~ is_low + is_high+ tt(is_high),
     data = bmt %>%
    mutate(is_low = ifelse(disgroup=="Acute Myelocytic Leukemia (Low Risk)",
                           1, 0),
           is_high = ifelse(disgroup=="Acute Myelocytic Leukemia (High Risk)",
                            1, 0)),
      tt=function(x,t,...) x*t )

tbl_regression(disgroup_mod, exponentiate=TRUE, conf.level=.9,
               label=list(is_low="Acute Myelocytic Leukemia (Low Risk)",
                          is_high="Acute Myelocytic Leukemia (High Risk)",
                          `tt(is_high)`="tt(Acute Myelocytic Leukemia (High Risk))"))

#-------------------------------------------------------------
# add time-varying covariates based on residual plots above 
#-------------------------------------------------------------

# piecewise linear interaction with time (to handle U-shapes)
# before the cutpoint, slope is beta_1
# after the cutpoint, slope is beta_2
# continuous at cutpoint 
# model_list[["mtx"]] <- coxph(
#   Surv(tdfs, deltadfs) ~ mtx + tt(mtx),
#      data = bmt %>% mutate(mtx = ifelse(
#        mtx=="Prophylactic Use of Methotrexate", 1, 0)),
#       tt=function(x,t,...) cbind( x * pmin(t, 600), x * pmax(0, t - 600))) 


model_list[["mtx"]] <- coxph(
  Surv(tdfs, deltadfs) ~ mtx + tt(mtx),
     data = bmt %>% mutate(mtx = ifelse(
       mtx=="Prophylactic Use of Methotrexate", 1, 0)),
      tt=function(x,t,...) x*t)

model_list[["disgroup"]] <- coxph(
  Surv(tdfs, deltadfs) ~ is_low + is_high+ tt(is_high),
     data = bmt %>%
    mutate(is_low = ifelse(disgroup=="Acute Myelocytic Leukemia (Low Risk)",
                           1, 0),
           is_high = ifelse(disgroup=="Acute Myelocytic Leukemia (High Risk)",
                            1, 0)),
      tt=function(x,t,...) x*t )

#---------------------------------------------------------------------------
# create combined table with HR's for each univariable model 
#---------------------------------------------------------------------------
model_summary_list <- pmap(list(var=covariates,
                                nice_name=nice_names, 
                                model=model_list),
                           get_model_summary) 

tbl_stack(model_summary_list, 
          group_header=paste0('Model: ', nice_names)) %>%
  as_gt() %>%
  gt::tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_row_groups()
  )

rm(model_summary_list)

# Fit model with continuous covariate
fit <- coxph(Surv(tdfs, deltadfs) ~ donorage, data = bmt)

# Choose representative values (e.g., quartiles)
quantiles <- quantile(bmt$donorage, probs = c(.1, .5, .9))
newdat <- data.frame(donorage = quantiles)

# Generate survival predictions
sf <- survfit(fit, newdata = newdat)

sf$surv %>%
  as_tibble() %>%
  mutate(time = sf$time) %>%
  pivot_longer(-c(time)) %>%
  ggplot(aes(x=time, y=value, color=name)) +
  geom_step()



tibble(time=sf$time,
       surv = sf$surv) %>%
  ggplot(aes(x=time, y=surv)) +
  geom_step()

plot(sf)

ggsurvfit(sf)

# perfect collinearity 
table(bmt$hospital, bmt$mtx)

bmt %>%
  filter(ts < ta) %>%
  nrow()

cox.zph(model, transform="km")$table %>%
  as.data.frame() %>%
  kbl(caption="Proportional Hazards Test") %>%
  kable_paper()

mod_resid <- as.data.frame(cox.zph(model, transform="km")$y) %>%
  as_tibble()

time <- cox.zph(model, transform="km")$time

df <- mod_resid %>%
  bind_cols(time=time)

df %>%
  pivot_longer(-time) %>%
  ggplot(aes(x=time,y=value)) +
  geom_point() +
  geom_smooth(formula=y~x,method="loess",span=5) +
  facet_wrap(~name, scales="free_y") +
  theme_c()

mod_resid <- as.data.frame(cox.zph(model, transform="km")$y) %>%
  as_tibble()

time <- cox.zph(model, transform="log")$time

df <- mod_resid %>%
  bind_cols(time=time)

df %>%
  pivot_longer(-time) %>%
  ggplot(aes(x=time,y=value)) +
  geom_point() +
  geom_smooth(formula=y~x,method="loess",span=5) +
  facet_wrap(~name, scales="free_y") +
  theme_c()


df %>%
  pivot_longer(-time) %>%
  ggplot(aes(x=time,y=value)) +
  geom_point() +
  geom_smooth(formula=y~x,method="loess",span=5) +
  facet_wrap(~name, scales="free_y") +
  theme_c()

model <- coxph(Surv(ta, deltaa) ~
                 mtx + age + donorage +
                 tt(waittime) + 
                 fab + disgroup+
                 sex_mismatch  +
                 cmv_mismatch, data=bmt,
               tt=function(x,t,...) x*splines::ns(t,df=2))

#------------------------------------------------------------
# g-computation to compute marginalized survival curves 
#------------------------------------------------------------
surv_treat <- survfit(model, newdata=bmt %>% 
          mutate(mtx="Prophylactic Use of Methotrexate")) %>%
  broom::tidy() %>%
  rowwise() %>%
  mutate(
    mean_surv = mean(c_across(contains("estimate")), na.rm = TRUE)
  ) %>%
  ungroup() %>%
  select(time, mean_surv) 

surv_no_treat <- survfit(model, newdata=bmt %>% 
          mutate(mtx="No Prophylactic Use of Methotrexate")) %>%
  broom::tidy() %>%
  rowwise() %>%
  mutate(
    mean_surv = mean(c_across(contains("estimate")), na.rm = TRUE)
  ) %>%
  ungroup() %>%
  select(time, mean_surv) 

surv_est_dir6 <- surv_treat %>%
  mutate(treatment="Prophylactic Use of Methotrexate") %>%
  bind_rows(surv_no_treat %>% 
              mutate(treatment="No Prophylactic Use of Methotrexate"))

surv_est_dir6 %>%
  ggplot(aes(x=time, y=mean_surv, color=treatment))+
  geom_step() +
  coord_cartesian(xlim=c(0,150)) +
  theme_c(axis.title=element_text(size=11),
          legend.text=element_text(size=11),
          legend.title=element_text(size=13)) +
  viridis::scale_color_viridis(
      discrete=TRUE, option="mako", end=.8, begin=.2) +
  labs(x="Time (Days)", y = "Event-Free Survival")

ggsave(here("figures/directive_6.jpeg"))

# cox.zph(model)


broom::tidy(model, 
            exponentiate=TRUE, 
            conf.int=TRUE) %>%
  mutate(term=gsub(
    "sex_mismatch|disgroup|cmv_mismatch|fab|mtx",
    "",term)) %>%
  mutate(term=factor(term),
         term=fct_reorder(term,estimate)) %>%
  ggplot(aes(y=term, xmin=conf.low,xmax=conf.high))+
  geom_point(aes(x=estimate), alpha=.4, size=1.1) +
  geom_errorbar(width=.2) +
  geom_vline(aes(xintercept=1),
             color='darkred', 
             linetype=2) +
  theme_c(axis.text.y=element_text(size=12, hjust=1)) +
  labs(y="", x= "Hazard Ratio")

gtsummary::tbl_regression(model,conf.level=.9,exponentiate=TRUE,
                           pvalue_fun = label_style_pvalue(digits = 3))

# ggsave(here("directive_6.jpeg"))

#--------------------------------------------------------
# obtain bootstrapped differences in survival curves
#--------------------------------------------------------

get_marginal_surv <- function(dat,
                              input_times = c(14, 30, 60, 90)) {
  
  model <- coxph(Surv(ta, deltaa) ~
                 mtx + age + donorage +
                 waittime + 
                 fab + disgroup+
                 sex_mismatch  +
                 cmv_mismatch, data=dat
               )

  surv_treat <- survfit(model, newdata=dat %>% 
            mutate(mtx="Prophylactic Use of Methotrexate")) 
  
  surv_treat_df <- tibble(time = input_times) %>%
    bind_cols(summary(surv_treat, times=input_times)$surv) %>%
    rowwise() %>%
    mutate(
      mean_surv = mean(c_across(!matches("time")), na.rm = TRUE)
    )  %>%
    select(time,surv_treat=mean_surv)
  
  surv_no_treat <- survfit(model, newdata=dat %>% 
            mutate(mtx="No Prophylactic Use of Methotrexate")) 
  
  surv_no_treat_df <- tibble(time = input_times) %>%
    bind_cols(summary(surv_no_treat, times=input_times)$surv) %>%
    rowwise() %>%
    mutate(
      mean_surv = mean(c_across(!matches("time")), na.rm = TRUE)
    )  %>%
    select(time,surv_no_treat=mean_surv)
  
  surv_no_treat_df %>% 
    left_join(surv_treat_df) %>%
    mutate(difference=surv_treat-surv_no_treat)
  
}
  

boot_surv <- map_df(1:5000, ~ {
  ind <- sample(1:nrow(bmt), replace=TRUE)
  boot_samp <- bmt[ind,]
  get_marginal_surv(boot_samp)
}) 

saveRDS(boot_surv, here('data/boot_surv_directive_6.RDS'))

est <- surv_est_dir6 %>% 
  pivot_wider(names_from=treatment, values_from=mean_surv) %>%
  mutate(difference=`Prophylactic Use of Methotrexate`-`No Prophylactic Use of Methotrexate`) %>%
  filter(time == 10 | time == 30 | time == 53 |
           time== 88) %>%
  pull(difference)


boot_surv <- readRDS(here('data/boot_surv_directive_6.RDS'))

boot_surv %>% 
  group_by(time) %>%
  summarize(lower = quantile(difference,.05),
            upper=quantile(difference,.95)) %>%
  bind_cols(estimate=est) %>%
  mutate(
         time_lab = paste0(time, " Days"),
         time_lab=factor(time_lab),
         time_lab=fct_reorder(time_lab,time)) %>%
  ggplot(aes(y=time_lab, xmin=lower,xmax=upper)) +
  geom_point(aes(y=time_lab,x=estimate))+
  geom_errorbar(width=.3) +
  theme_c(axis.title.x=element_text(size=16),
          axis.text.y=element_text(size=14)) +
  geom_vline(aes(xintercept=0), color="darkred", linetype=2) +
  labs(x = latex2exp::TeX("$\\hat{S}_1(t) - \\hat{S}_0(t)$"),
       y="") 

library(kableExtra)

boot_surv %>% 
  group_by(time) %>%
  summarize(lower = quantile(difference,.05),
            upper=quantile(difference,.95)) %>%
  mutate(across(where(is.numeric), ~round(.x,4))) %>%
  kbl(caption="Bootstrap Confidence Intervals for $\\hat{S}_1(t) - \\hat{S}_0(t)$") %>%
  kable_paper()

# Age, fab classification, wait time until transplant, donor age, donor sex, cmv, sex, donor cmv, disease group

dat <- bmt %>% 
  select(age,fab, waittime, disgroup, 
         deltap, tp, deltadfs, sex_mismatch,
         cmv_mismatch,donorage,mtx,
         tdfs) %>%
  mutate(id=row_number())

dat <- tmerge(dat, dat, id=id,endpt=event(tdfs,deltadfs))
dat <- tmerge(dat, dat, id=id, plat=tdc(tp, deltap,init=0))

model <- coxph( Surv(tstart, tstop, endpt) ~ age + 
                  #waittime + fab +
         cmv_mismatch  + disgroup + plat +donorage + sex_mismatch +
           mtx, data=dat) 

gtsummary::tbl_regression(model,
                          label=list(plat="Recovery of Normal Platelet Levels",
                                     age="Age",
                                     fab="FAB Classification",
                                     male="Sex",
                                     disgroup="Disease Group",
                                     cmv="Patient CMV",
                                     donorage="Donor Age",
                                     sex_mismatch="Sex Mismatch",
                                     mtx="Treatment" #,
                                     # waittime="Wait Time"
                                     ),
                          exponentiate=TRUE,
                          conf.level=.9,
                           pvalue_fun = label_style_pvalue(digits = 3))

dat <- bmt %>% 
  select(age,fab, waittime, disgroup, 
         deltap, tp, deltadfs, deltar,
         deltas, sex_mismatch,cmv_mismatch,donorage,
         mtx,
         tdfs) %>%
  mutate(id=row_number())

dat <- tmerge(dat, dat, id=id,endpt=event(tdfs,deltadfs))
dat <- tmerge(dat, dat, id=id, plat=tdc(tp, deltap,init=0))



library(timereg) 

dat <- dat %>%
  mutate(cause = case_when(
    deltas == 1 & deltar == 0 ~ 1,
    deltar == 1 ~ 2,
    deltar == 0 & deltas == 0 ~ 0
  ),
  cause=as.integer(cause))

bmt_copy <- bmt

odds_mod <- prop.odds.subdist(Event(tdfs, cause) ~ plat, 
         cause = 2,
         data=dat)

odds_mod <- prop.odds.subdist(Event(tdfs, cause) ~ 
                           age + # fab+ waittime+# fab +# waittime + 
         cmv_mismatch  + disgroup  +donorage + sex_mismatch +
           plat + mtx, 
         cause = 2,
         data=dat)

coef(odds_mod) %>% as.data.frame() %>%
  as_tibble(rownames="variable") %>%
  mutate(across(where(is.numeric), ~round(.x,4))) %>%
  kbl() %>%
  kable_paper()

coef(odds_mod) %>% as.data.frame() %>%
  as_tibble(rownames="variable") %>%
  rename(coef=`Coef.`) %>%
  mutate(lower = coef + qnorm(.05)*SE,
         upper= coef+ qnorm(.95)*SE) %>%
  mutate(across(c(coef,lower,upper),exp)) %>%
  mutate(variable=gsub(
    "sex_mismatch|disgroup|cmv_mismatch|fab|mtx",
    "",variable),
    variable=gsub("plat", "Recovery of Normal Platelet Levels", variable),
    variable=gsub("waittime", "Wait Time for Transplant", variable),
    variable=gsub("age", "Age", variable),
    variable=gsub("donorAge", "Donor Age", variable)) %>%
#  mutate(across(c(`Coef.`, `lower2.5%`,`upper97.5%`), exp)) %>%
  ggplot(aes(y=fct_reorder(variable,coef),
             xmin=lower,
             xmax=upper)) +
  geom_errorbar(width=.2) +
  geom_point(aes(x=coef),size=.55) +
  geom_vline(xintercept=1,linetype=2,color="darkred") +
  theme_c(axis.text.y=element_text(size=11, hjust=1),
          axis.title.x=element_text(size=11.5)) +
  labs(y="", x="Odds Ratio")

ggsave(here("figures/directive_7.jpeg"))



coef(odds_mod) %>% as.data.frame() %>%
  as_tibble(rownames="variable") %>%
  rename(coef=`Coef.`) %>%
  mutate(lower = coef + qnorm(.05)*SE,
         upper= coef+ qnorm(.95)*SE,
         lower=exp(lower),
         upper=exp(upper),
         coef=round(exp(coef),3),
         CI= paste0(round(lower,3), ", ", round(upper, 3)),
         `P-val`= round(`P-val`,3))  %>%
  select(Characteristic=variable, OR=coef, `90% CI`=CI,
         `p-value`=`P-val`) %>%
  mutate(Characteristic = gsub(
    "cmv_mismatch|disgroup|sex_mismatch", "", 
    Characteristic),
    Characteristic=str_to_title(Characteristic),
    Characteristic=gsub("Cmv", "CMV", Characteristic),
    Characteristic=gsub("Donorage", "Donor Age", Characteristic),
     Characteristic=gsub("Plat", 
                         "Recovery of Normal Platelet Levels",
                         Characteristic)) %>%
  kbl() %>%
  kable_styling(bootstrap_options="striped")